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We describe a Monte Carlo event generator for the simulation of QCD-instanton 
induced processes in deep-inelastic scattering (HERA). The QCDINS package is 
designed as an "add-on" hard process generator interfaced to the general hadronic 
event simulation package HERWIG. It incorporates the theoretically predicted pro- 
duction rate for instanton-induced events as well as the essential characteristics 
that have been derived theoretically for the partonic final state of instanton-induced 
processes: notably, the flavour democratic and isotropic production of the partonic 
final state, energy weight factors different for gluons and quarks, and a high average 
multiplicity O(10) of produced partons with a Poisson distribution of the gluon mul- 
tiplicity. While the subsequent perturbative evolution of the generated partons is 
always handled by the HERWIG package, the final hadronization step may option- 
ally be performed also by means of the general hadronic event simulation package 
JETSET. 
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PROGRAM SUMMARY 



Title of program: QCDINS 2.0 
Catalogue identifier: 



Program obtainable from: http:/ / www.desy.de/~tOOfri/qcdins/qcdins.htm] 

Computer for which the program is designed and others on which it has been tested: 
Any computer with a FORTRAN 77 compiler 

Operating systems under which the program has been tested: Linux 2.0.X; HP-UX 
10.2 



Programming language used: FORTRAN 77 

Memory required to execute with typical data: Size of executable program is approxi- 
mately 2.6 MB. The size of the QCDINS library itself is about 200 KB; the required 
routines from the HERWIG and JETSET libraries constitute the dominant portion 
of the needed memory. 

No. of processors used: 1 

Has the code been vectorised or parallelized?: no 

No. of bytes in distributed program, including test data, etc.: 1071106 
Distribution format: ASCII 

CPC Program Library subprograms used: HERWIG [1] version 5.9; JETSET 7.4 [2] 
Keywords: QCD; Instanton; Deep-inelastic scattering; Monte Carlo simulation 
Nature of physical problem 

Instantons are a basic aspect of Quantum Chromodynamics. Being non-perturbative 
fluctuations of the gauge fields, they induce hard processes absent in conventional 
perturbation theory. Deep-inelastic lepton-nucleon scattering at HERA offers a 
unique possible discovery window for such processes induced by QCD-instantons 
through their characteristic final-state signature and a sizable rate, calculable within 
instanton-perturbation theory. An experimental discovery of such a novel, non- 
perturbative manifestation of non-abelian gauge theories would be of fundamental 
significance. However, instanton-induced events are expected to make up only a 
small fraction of all deep-inelastic events. Therefore, a detailed knowledge of the re- 
sulting hadronic final state, along with a multi-observable analysis of experimental 
data by means of Monte Carlo techniques, is necessary. 



Method of solution 

The QCDINS package is designed as an "add-on" hard process generator interfaced 
to the general hadronic event simulation package HERWIG. It incorporates the 
theoretically predicted production rate for instanton-induced events as well as the 
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essential characteristics that have been derived theoretically for the partonic final 
state of instanton-induced processes: notably, the flavour democratic and isotropic 
production of the partonic final state, energy weight factors different for gluons 
and quarks, and a high average multiplicity O(10) of produced partons with a 
Poisson distribution of the gluon multiplicity. While the subsequent perturbative 
evolution of the generated partons is always handled by the HERWIG package, the 
final hadronization step may optionally be performed also by means of the general 
hadronic event simulation package JETSET. 

Restrictions on the complexity of the problem 

The default values of the implemented kinematical cuts represent the state of the 
art limits for the reliability of the generated instanton-induced event rate and event 
topology. 

Typical running time 

10 - 100 events per second for a PC with Pentium CPU, depending on its clock 
frequency. On a HP 9000/735 (99 MHz) workstation, 6 events per second are gen- 
erated. 

Unusual features of the program 
none 

References 

[1] G. Marchesini et al., Comput. Phys. Commun. 67 (1992) 465. 
[2] T. SjSstrand, Comput. Phys. Commun. 82 (1994) 74. 
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LONG WRITE-UP 



1 Introduction 

The ground state ( "vacuum" ) of non-abelian gauge theories like QCD is known 
to be very rich. It includes topologically non-trivial fluctuations of the gauge 
fields, carrying an integer topological charge Q. The simplest building blocks 
of topological structure in the vacuum are [1,2] instantons with Q — +1 and 
anti-instantons with Q — — 1. Instantons represent gluon field configurations 
that are localized ( "instantaneous" ) in Euclidean time and space. While they 
are believed to play an important role in various long-distance aspects of QCD, 
there are also important short-distance implications. In QCD with n/ (light) 
flavours, instantons induce hard processes violating chirality in accord [2] with 
the selection rule A chirality = 2nfQ, due to the general chiral anomaly. While 
in ordinary perturbative QCD (Q = 0) these processes are forbidden, their ex- 
perimental discovery would clearly be of basic significance. The deep-inelastic 
scattering regime is strongly favoured in this respect, since hard instanton- 
induced processes are both calculable [3,4] within instanton-pertubation the- 
ory and have good prospects for experimental detection at HERA [4-7]. 

QCDINS [6] is a Monte Carlo package for simulating QCD-instanton induced 
scattering processes in deep-inelastic scattering (HERA). It is designed as an 
"add-on" hard process generator interfaced by default to the Monte Carlo 
generator HERWIG [8]. It incorporates the theoretically predicted production 
rate for instanton-induced events as well as the essential characteristics that 
have been derived theoretically for the partonic final state of instanton-induced 
processes: notably, the flavour democratic [2] and isotropic [5,9] production of 
the final state partons, energy weight factors different for gluons and quarks [3], 
and a high average multiplicity 2rif + 0(l/a s ) of produced partons with a 
(approximate) Poisson distribution of the gluon multiplicity [5,9,10]. 

Earlier versions of QCDINS have been used already to establish first experi- 
mental bounds on the rate of instanton-induced events at HERA [11-13] and 
to develop instanton search strategies [7]. 

In the present report a comprehensive description is given of the theoretical 
framework built into the program (Section 2) as well as of the various program 
components (Section 3) and of their usage (Section 4). 
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2 Instanton-induced events in deep-inelastic scattering 



Let us briefly summarize in this section the underlying physics picture, some 
relevant formulae and the main stages involved in QCDINS to generate the 
complete instanton-induced partonic final state. The remaining formulae may 
be found under the corresponding descriptions of QCDINS routines in Sec- 
tion 3. 



e 




Fig. 1. Structure and kinematic variables of the dominant instanton-induced process 
in deep-inelastic scattering. 

In deep-inelastic e ± P scattering, instanton-induced events are predominantly 
associated [4] with a process structure as sketched in Fig. 1: A photon splitting 
into a gg-pair fuses with a gluon from the proton in the background of an in- 
stanton (/) or an anti-instanton (/). For each (light) flavour, q — {d, u, s, . . .}, 
a violation of chirality is induced, 



A chirality (q) = A [#(q R + q R ) - #(<?l + q L )\ = ± 2 for an < 



_ , (1) 
I 



in agreement with the general chiral anomaly relation [2]. Correspondingly, the 
partonic final state exhibits "flavour democracy", i.e. qRq R (qlQl) P a hs of all 
rif light flavours occur precisely once in case of an instanton (anti-instanton), 
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As illustrated in Fig. 1, one of those partons acts as a current-quark (jet) k, 
whereas the other 2n/ — 1 (anti-)quarks and some number n g of gluons are 
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directly emitted from the instanton (anti-instanton) "blob" . Instanton-induced 
processes initiated by a quark from the proton are suppressed by a factor of 
a 2 s with respect to the gluon initiated process [4]. This fact, together with the 
high gluon density in the relevant kinematical domain at HERA, justifies to 
neglect quark initiated processes. 

In instanton-perturbation theory, the dominant instanton-induced contribu- 
tion to the inclusive eP cross sectionF], subject to appropriate kinematical 
restrictions and (theoretical) fiducial cuts, has a convolution-like structure [4], 
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It involves integrations over the gluon density f g (z), the virtual photon flux 
and the flux of virtual (anti-)quarks q' in the instanton-background [4,14], 

^ » £SL(i + _L_J._ ( 4 ) 

16 7T J ZX' \ X B j X' bXBjyBjJ 



All relevant kinematical variables in Eq. (3) are defined in Fig. 1, and A 
denotes the electric charge squared of the virtual (anti-)quark q' in units of 
the electric charge squared e 2 = Ana. 

In Eq. (3), a^) denotes the instanton-induced total cross section of the q'g- 
subprocess (c. f. Fig. 1) and contains the essential instanton dynamics. Its 
analytical form [4] used in QCDINS may be found in Eq. (17). As illustrated 

3 A sum over instanton-induced and anti-instanton induced processes is always 
implied by the superscript (I) at cross sections. 
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Fig. 2. Instanton-subprocess cross section (17) from Ref. [4] for n/ = 3 and A^ 
from Eq. (7). It exhibits a strong dependence on the corresponding Bjorken variables 
Q' 2 and x', respectively. 

in Fig. 2, a^? g is very steeply growing for decreasing values of Q' 2 and x', respec- 
tively. Eventually, our theoretical predictions based on instanton-perturbation 
theory [4] will cease to hold. Therefore, the following cuts inferred from a high- 
quality lattice simulation of QCD [15] are implemented by default in QCDINS 
2.0 (Table 4), 



g / >g; Ilin = 3o.8A 



MS ' 



x <&Yi\\Yi — 0.35. 



(5) 



A further cut on the photon virtuality, 



Q = JSxBjVBi > Qmin = Q' m 



(6) 



is applied in order to warrant sufficient suppression of non-planar contribu- 
tions [3], which are hard to calculate and may spoil the validity of Eq. (3). 

The cross section a^ g exhibits a rather weak residual dependence on the renor- 
malization scale /i r . As an "optimal" choice, the value fi r = 0.15 Q', corre- 
sponding to the minimum [4], da q i g /dfi r = 0, is taken by default (Table 4). 



However, a^Sl depends strongly on the QCD scale A.^' . Since, strictly speak- 



er 



MS 



ing, the underlying theoretical framework refers to massless quarks, the (de- 
fault) number of flavours is set to rif = 3 (Table 4). The respective value 



Ag> = 0.346 +0 -° 31 GeV, 
MS -0.029 



(7) 



is obtained by a standard 3-loop perturbative flavour reduction (Eq. (9.7) of 
Ref. [16]) from the 1998 world-average of the running QCD coupling at the 
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Z-boson mass [16], 



a s (M z ) = 0.119 ± 0.002 



3 1T A§ = 0.219 



+0.025 
-0.023 



GeV. 



(8) 



The central values of these parameters are taken as default in QCDINS 2.0 
(Table 4). This upgrade of together with the modified cuts (5) and (6), 

represents an improved understanding of the input parameters and a consid- 
erable reduction of uncertainties, as compared to the original publication [4] 
and earlier versions of QCDINS. Note that it also implies a significant change 
in the predicted rate. 

Next, let us summarize the various stages of event generation by means of 



In a first stage, the various Bjorken variables Q' 2 , x', z, XBj, |/Bj of the instanton- 
induced process (c. f. Fig. 1) are generated, with a distribution according to 
the normalized differential cross section from Eq. (3), 



In the second stage of the event generation, the 4- momenta g, q, q', k of the 
incoming gluon g, the virtual photon q, the virtual quark q' and the current 
quark k, respectively, are filled. Sudakov decompositions of these momenta 
are used to incorporate various constraints, e.g. on the momenta squared. The 
4-momentum e' of the outgoing lepton is calculated subsequently. 

In the third stage, the partonic final state of the instanton-induced g'g-subpro- 
cess is generated in its centre-of-mass system (CMS) as follows. The number 
n g of produced gluons is generated according to a Poisson distribution with 
mean {n g )^\x' ,Q') ~ l/a s ~ 3 (for the cuts (5)), as calculated theoretically 
(Eq. (13)) in instanton-perturbation theory [5,9,10]. Next, n/(= 3) [q . . . q]- 
"strings" of partons are set up, each beginning with a quark, followed by a 
random number < n g + 1 of gluons and ending with an anti-quark of randomly 
chosen flavour. There are n g + 1 gluons in total and, due to the required flavour 
democracy (2), gg-pairs of all n/ flavours occur precisely once. A quark and a 
gluon among these 2n/ + n g + 1 partons are (randomly) marked as incoming. 

The momenta Pi of the n = 2rif — 1 + n g outgoing partons are then generated 
in the CMS of the instanton subprocess, according to the energy-weighted 
phase-space 



QCDINS. 



1 d b <j { J), 



(9) 



^ dQ' 2 dx' dz dxBj efa/Bj 
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These different energy weights for quarks and gluons [3], along with the angular 
isotropy [5,9], are characteristic features of the leading-order partonic final 
state (after averaging over colour). 

Next, the colour and flavour connections of the partons are set up. The colour 
flow is obtained simply by connecting the colour lines of adjacent partons 
within each of the above-mentioned n/ [q . . . q]- "strings" in a planar manner 
(consistent with the leading order 1/N C expectation). This choice is inspired by 
the leading-order partonic final state (after averaging over colour) [3] , but may 
well deserve further research. The flavour flow is constructed by connecting 
the flavour lines of the quark at the beginning of a string with the flavour line 
of the anti-quark at the end of a string. 

The hard subprocess generation ends by boosting the momenta of the final 
state partons to the laboratory frame. 

While the subsequent perturbative evolution of the generated partons is al- 
ways handled by the HERWIG [8] package, the final hadronization step may 
optionally^ be performed also by means of JETSET [18]. 



3 The QCDINS package 

This section is devoted to a systematic description of the various routines of 
the QCDINS package that is designed as an "add-on" hard process generator, 
interfaced to the Monte Carlo generator HERWIG [8]. 

This reference section is organized as follows: 

While all subroutines and functions of the QCDINS package are described in 
alphabetical order in Section 3.1, a logical flow-chart is provided in form of 
Tables 1 and 2 below. They should always be used as the main guide through 
the description of the package. A routine listed in the n-th column and the 
m-th row of these tables progressively calls all routines in the (n+l)-th column 
starting in the (m+l)-th row. All routines called at the level of the main hard 
process generator QIHGEN and below are documented in Table 2. 

4 We thank H. Jung for his help [17] to interface QCDINS with JETSET. 
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Initialize input para- 
meters (Tables 4, 5) 


QIINIT 










Flavour reduction: 

A (5) =>A W< 5 ) 
MS MR 




LAMNF 








±111 l)ld'llZ/\_. l_/CXl UlL^lC 

INST in event record 




QIINIH 








P vitri" mmif Q V VY"1 ^ 
JT 1 111JJU.L ^Jctl Cl.Illt. 

-j-pT-q Q ti fl w£l VT1 1 Tl 0"<> 

Lr~l J (XllVA VV Cl± 11111 tiO 


QISTAT 










Initialize JETSET 
common block data 


GJEINI 










Loop over desired 
TininhpT of" pvpnt^i 

llLllllk_/^± Wl 


QCLOOP 










Generate one event 




QCDGEN 








A eel (r r i norn nrru^r>cc 

JT-SSlgll Ildl U piULLM 

variables (HERWIG) 






HWEPRO 






No action (modified 
HERWIG routine) 








HWEGAM 




Call main instanton 
process generator 








HVHBVI 




Main (hard) instanton 
process generator 










QIHGEN 


frOTir'T'nt'f 1 Tin f(~(~iTi 

VI L 1 1<_ 1 (X L JJ<X1 LrWll 

cascades (HERWIG) 






HWBGEN 






Combine jets with 
correct kinematics 








HWBJCO 




Convert HERWIG to 
JETSET block data 






HERLUND 






JETSET event record 
to HEPEVT common 






LUHEPC 







Table 1 

Flow-chart of QCDINS routines: A routine in the n-th column and the m-th row 
progressively calls all routines in the (n+l)-th column starting in the (m+l)-th row. 
All routines called at the level of the main hard process generator QIHGEN and 
below are listed in Table 2. 

A specific application requires the writing of a steering program by the user 
(c. f. Section 4 and Appendix A). It must contain the standard HERWIG- 
initialization calls as well as the calls to various initialization routines for 
QCDINS. The latter comprise essentially the four routines listed in the second 
column of Table 1. By calling the last one of these (QCLOOP), the user starts 
the proper simulation which comprises the chain of internally called QCDINS 
routines as documented in Tables 1 and 2. A specific and quite extensive 
example is provided with the QCDINS distribution ('qtesthz.F') and may be 
found in the directory 'qedtest'. It also illustrates the use of the event analysis 
routine HWANAL called by HERWIG after each processed event. 
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Optional account of initial state 
radiation from lcpton 


EXFRAC 






Check kinematical boundaries 


QICALC 






Generate identity code of current 
quark k and virtual quark q' 


QIHPAR 






Generate (Q',x') and associated weight 


QIHINS 






Generate X = Q' 2 ,x' as dX/X^+V 




QIRDIS 




Calculate total cross section of 
instanton-induccd subprocess q'g =>' X 




Q2SIG 




Calculate //-valley action 






ACTION 


Calculate termionic overlap ui(£) 






OMEGA 


Calculate Lambert W-function 








Calculate saddle-point value of 
conformally invariant //-distance £ 






XI 


Calculate inverse running coupling 4 "~ 

"ms 






XQS 


Generate number n ff of emitted gluons 








Calculate gluon multiplicity (n g )^ 




1\ /TT TT r l'' 




Calculate inverse running coupling 4 "~ 

"ms 






XQS 


Calculate /7-valley action S {II ~>(0 






ACTION 


Calculate Lambert W-function 






LAMBERT W 


Calculate flux of virtual quark q' 


QISPLT 






Calculate remaining weight 


QIPVWT 






Calculate momentum of incoming gluon 


QIKPAR 






Generate momentum of virtual photon 


QIKGAM 






Generate momenta of virtual quark q' 
and current quark k 


QIKGSP 






Generate partonic final state 


QISTID 






Generate the partons of the instanton 
subprocess in form of [q . . . q]-" strings" 




QIGLST 




Find the incoming partons 
in the [q . . . g]-strings 




QIGPAR 




Assign masses of outgoing partons 




QIPLST 




Generate 4-momenta of outgoing partons 




QIPSGN 




Calculate relative energy weight 
of outgoing partons 






QIPSWT 


Store 4-momenta of outgoing partons 
into PHEP common block of HERWIG 




QIPSTO 




Colour/flavour connections for each string 




QICCON 





Table 2 

All routines called at the level of the main hard process generator QIHGEN and 
below. As in Table 1, a routine in the n-th column and the m-th row progressively 
calls all routines in the (n+l)-th column starting in the (m+l)-th row. The displayed 
sequence of calls corresponds to the default settings of control flags. 
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3. 1 Subroutines and functions 



SUBROUTINE ACTION (XI4, S, DS, DDS) 

Purpose: Calculation of the //-action as well as its 1st and 2nd derivatives, 
as function of the conformally invariant //-distance. 



Arguments: 

XI4: conformally invariant //-distance £. 

S: //-action 5^(0, Eq. (11). 

DS: dS^/dt, 

DDS: d 2 S^/de 

Procedure: 



The action is calculated according to the exact valley form [19,20], 



5^(0 = 1- + 8] In 



/(0 /(0 2 /(0 3 
/(0 = e 2 + v / ? T 4e- 4 - 



^(/(0 + 4) 



(11) 
(12) 



SUBROUTINE EXFRAC (A) 

Purpose: Optional account of initial state radiation from the lepton. 

Arguments: 

A: = 1.0; dummy rescaling factor of the incoming lepton mo- 
mentum; double precision output variable. 

Remarks: The actual routine has to be provided by the user. 
SUBROUTINE GJEINI 

Purpose: Initialization of the JETSET [18] parameter common blocks LU- 
DAT1, LUDAT2, LUDAT3, LUDAT4 and LUDATR. 



Remarks: GJEINI has to be called by the user before any other JETSET 
routine. GJEINI is from Ref. [21]. 
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FUNCTION 

GMULT (XPR, XLMIN, XLMAX, QLAM, KAPPA, NF, LOOPFL) 



Purpose: 



Calculation of the average gluon multiplicity (n g )^ depending 

on x', Q'/A^\ iir/Q'i n f and loop-order. Here, p r and rif de- 
note the renormalization scale and the number of light flavours, 
respectively. 



Arguments: 



Procedure: 



XPRIME: x' 
XLMIN: £mi 
XLMAX: £ ma 
QLAM: Q'/A^' : 
KAPPA: 
NF: 



i, lower boundary of £ used for interpolation. 
x ; upper boundary of £ used for interpolation. 



MS 



Hr/Q' 

rif] number of light flavours. 



LOOPFL: =1: 1-loop renormalization group (RG) invariance [4] 
along with 1-loop form of a s . 

= 2: 2-loop RG invariance [4] along with 2-loop form 



of a m . 



3: (default) 2-loop RG invariance along with 3-loop 



form of a^g. 



From an analysis based on the generalized (Mueller [22]) optical theorem for 
the q'gg forward scattering amplitude and the //-valley method, one infers [9] 
the differential one-gluon inclusive q'g 4> g + X cross section, normalized by 



the total cross section a^ g . The mean gluon multiplicity [5,10] is then found 
by phase space integration, 



CO 



x 



Q' Hr 



2tt 



' A^£ } ' Q' 

MS 



n f = 



CA- 



MS 



(1/p*) 



c, ^dS™,^ 

(£*- 2 )^HU 



(13) 



The function GMULT calculates and returns the average gluon multiplic- 
ity (13). 

The stars (*) in Eq. (13) denote the saddle point values of the // collective 
coordinates p, p and £. Their computation proceeds as in the descriptions of 
the functions Q2SIG and XI. The required values of 47r/o;j^g are calculated 
and returned by the function XQS. The //-action S( ir )(£) and its ^-derivatives 
are provided by the subroutine ACTION. 
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SUBROUTINE HERLUND 
Purpose: Conversion of the HERWIG 



Remarks: 



event record in the HEPEVT 
common block to the respective JETSET [18] common block. 

HERLUND is a modified [17] copy from the JETSET subrou- 
tine LUHEPC. 



SUBROUTINE HVHBVI 

Purpose: Call of the main (hard) instanton process generator QIHGEN. 

Remarks: Replaces a dummy stub in HERWIG [8] that was originally used 
as event generation interface for the Monte Carlo generator HER- 
BVI [23] for baryon number violating interactions. Used in the 
QCDINS package to select QCD-instanton induced processes via 
MOD (IPROC/100,100) > 75. The process code IPROC (= 17600) 
has to be set in the user's steering program (c. f. Appendix A). 
Furthermore, the QCDINS program header is printed. 

SUBROUTINE HWBJCO 

Purpose: Modification of HERWIG [8] routine to account for instanton- 
induced scattering. 

Remarks: The modifications are [24]: The logical flag (DISPRO) for keeping 
the lepton momenta fixed in HERWIG 5.9 is modified to include 
also instanton-induced DIS, IPRO = 76. Furthermore, a bug con- 
cerning energy-momentum conservation in the original routine 
from HERWIG 5.9 has been fixed. 

SUBROUTINE HWEGAM (IHEP, ZMI, ZMA, WWA) 

Purpose: Modification to avoid standard generation of the (virtual) photon 
at this stage. 

Arguments: c. f. Ref. [8] 

Remarks: This is a modified routine from HERWIG 5.9 (c.f. Ref. [8]). Usu- 
ally, HWEGAM generates an incoming photon from the incom- 
ing e^. Within QCDINS, however, the photon is generated at 
a later stage in the subroutine QIKGAM. Thus, the HERWIG 
routine HWEGAM has been modified to immediately return for 
instanton-induced processes (IPRO=76). 

FUNCTION LAMBERTW (X) 

Purpose: Calculation of the principal branch of the Lambert W-function 
W(x) for x > 0. 
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Arguments: X: x > 0; argument of the Lambert W-function. 
Procedure: 

W(x) is the (real) solution of W(x) exp(W(x)) = x, analytic at x — 0. The 
following simple, but accurate approximation is used and returned by LAM- 
BERTW: 



W(x) 



0.665 • (1 + 0.0195 • ln(x + 1)) ln(x + 1) + 0.04; 

for < x < 500; (14) 

k \n{x - 4) - (1 - j^y) • ln(ln(x)); for x > 500 . 



FUNCTION LAMNF (NF, LAMBDA5) 

Purpose: Calculation of A^- <5 ' ) from to order a^. 

r MS MS MS 

Arguments: 

NF: number of (light) flavours, rif. 

LAMBDA5: input value A^. 

Procedure: The flavour reduction of A^ to the desired number of light 
flavours is performed by using Eq. (9.7) of Ref. [16]. 

SUBROUTINE LUHEPC (MCONV) 

Purpose: Conversion of the JETSET [18] event record contents back to the 
HEPEVT common block. 



Arguments: MCONV = 1 

Remarks: The present routine is a modified [21] version of the JETSET 
routine LUHEPC. 

FUNCTION OMEGA (XI4) 

Purpose: Calculation of the fermionic overlap, as function of the confor- 
mally invariant //-distance. 



Arguments: XI4: conformally invariant //-distance £. 

Procedure: 



The following simple, but accurate approximation for the fermionic overlap [4] 
uj is used and returned by OMEGA: 

"(0 - (e + i 4 / 2 )3/ 2 - (15) 
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FUNCTION Q2SIG (XPRIME, QLAM, KAPPA, LOOPFL, NF) 
Purpose: Calculation of the total cross section 



Q' 2 a { Jl [nb GeV 2 ] (16) 

for the instanton-induced subprocess, depending on x', Q'/A 1 ^ , 
Hr/Q'i loop-order and rif. Here, \x r and rif denote the renormal- 
ization scale and the number of light flavours, respectively. 

Arguments: 

XPRIME: x' 
QLAM: Q'/A^ 
KAPPA: n r /Q' 

LOOPFL: =1: 1-loop renormalization group (RG) invariance [4] 
along with 1-loop form of a s . 
= 2: 2-loop RG invariance [4] along with 2- loop form 
of a m . 

= 3: (default) 2-loop RG invariance along with 3-loop 
form of a^g. 
NF: rif; number of light flavours. 

Procedure: 

The function Q2SIG calculates and returns the cross section (16) as derived 
in Ref. [4], 



n ,2 (/) _ ,2 v 7 ^ 15/2 ((£* + 2K 2 + 4S(S-2v.)) ( &-2) A^A 



7/2 



X ■ 



3„,* 5 



\{S - u„ - 2D(S)y + S(S - v*)D ( In (-^gL 



(17) 



x 



4tt 



19/2 



cxp 



47T 



"Ms M 



(&) - 2 1 - In 



It is expressed entirely in terms of the saddle point values for the II collective 
coordinates, £* (conformally invariant distance) and v* = Q'p* (scaled size). 
For given x', Q' / A^^ and (scaled) renormalization scale fi r /Q', these are in 
turn unique solutions of the saddle point equations [4] 



£*-2 = 4 



x 

l-x' 



1 - 



(18) 
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v. = 2D(S(&))W 



4tt 



1 , gtf*) 



V 



2I>(5(^)) 



(19) 



with 



and 



Ai 



1 for LOOPFL = 1; 

1 + for LOOPFL = 2, 3; 



(20) 



Ao = ^ 



12 ( 3o 2Mdtd for LOOPFL 



for LOOPFL = 1; 

= 2,3. 



(21) 



The //-action S( IT \£) as well its ^-derivatives, entering the cross section (17) 
and Eqs. (18), (19) through 

S&) = A 1( 3 SW(Q-A 2 ; £>(/(&)) = —jL—ffo) , (22) 



are calculated in the subroutine ACTION. The fermionic overlap u;(£) is 
calculated and returned by the function OMEGA. In Eq. (19), W denotes 
the principal branch of the Lambert lU-function, i.e. the (real) solution of 
W(x) exp(W(x)) = x, analytic at x — 0. The latter is calculated and returned 
by the function LAMBERTW. 

The first step in the solution of the saddle-point equations (18), (19) consists in 
eliminating v* in Eq. (18) by inserting Eq. (19). Next, for given x', Q' / A^^ and 
Hr/Q'i t ne resulting implicit equation is solved numerically for £*. This is done 
by the function XI which provides = Q'/A^-t , fj, r /Q') on return. The 

latter is then inserted into Eq. (19) providing i>* = v*(x', Q' /A^\ fi r /Q'). 

The values of 4ir/ay^ are calculated and returned by the function XQS. 

SUBROUTINE QCDGEN 

Purpose: Interface for the generation of one instanton-induced event, in- 
cluding calls to event initialization, hadronization and event ter- 
mination routines. 



Procedure: 
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• An instanton-induced event is initialized by a call to the HERWIG [8] sub- 
routine HWUINE. 

• The partonic instanton subprocess is generated by a call to the HERWIG 
subroutine HWEPRO. 

• If the hadronization flag QICONT(21) is set .TRUE, (default), the event is 
fully hadronized. Else, the event is finalized (HWUFNE) immediately after 
the call of HWEPRO. Depending on the control flag for the hadroniza- 
tion model, QICONT(18)=.TRUE./.FALSE., the hadronization is either 
performed by appropriate HERWIG [8] or (modified [17]) JETSET [18] 
routines, respectively (see Table 3). 

• Furthermore, in this routine energy-momentum conservation in the gener- 
ated event is checked. 

SUBROUTINE QCLOOP 

Purpose: Loop over calls to the instanton event generator and interface of 
the instanton-induced hard subprocess to HERWIG [8] initializa- 
tion and output routines. 



Procedure: 



• The instanton-induced hard subprocess is initialized by calling the HERWIG 
subroutine HWEINI. 

• The generator QCDGEN for an instanton-induced event and, thereafter, 
the user's routine HWANAL for analyzing data from the event, are called 
in a loop MAXEV times, with MAXEV being the desired number of events. 
MAXEV has to be set in the user's steering program (c. f. Appendix A). 

• The instanton-induced hard subprocess is terminated by calling the HER- 
WIG subroutine HWEFIN that produces the final output. 

SUBROUTINE QICALC 

Purpose: Calculates derived parameters and checks the chosen kinematical 
boundaries. The latter are adjusted in case of inconsistencies. 



Non-local variables (re)set: 
QILPAR(I): 



QIDPAR(l): 
QIWARF(l): 



= LOG(QIUPAR(I)), 1 = 1.. 18, if QIUPAR(I) 
> 1.0 • 1(T 20 ; otherwise QILPAR(I) = -1.0 • 
10 25 ; these parameter settings are partially re- 
set after exit (see below). 

= QIUPAR(12)**2; square of minimum allowed 
instanton CM energy Wf min . 
= .FALSE, (default); if <2 m ax violates the kine- 
matical constraint 



Q 



max — -H)ax* t) ' 



W? ■ 

" " 1 mm ' 



(23) 
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yiUUlN 1 (10 ) 


Hadronization scheme 




HERWIG Interface 




HWUINE 


Initialize event 




HWEPRO 


Generate hard subprocess 




HWBGEN 


Generate parton cascades 




HWDHQK 


Do heavy quark decays 


.TRUE. 


HWCFOR 


Do cluster hadronization 


HWCDEC 


Do cluster decay 




HWDHAD 


Do unstable particle decays 




HWDHVY 


Do heavy flavour decays 




HWMEVT 


Add soft underlying event if needed 




HWUFNE 


Finalize event 




JETSET Interface 




HWUINE 


Initialize event 




HWEPRO 


Generate hard subprocess 




HWBGEN 


Generate parton cascades 




HWUFNE 


Finalize event 




HERLUND 


Convert HERWIG event record in HEPEVT 
common block to the respective JETSET 
common block 


.FALSE. 


LUGIVE 


Assign the hard subprocess (HARD) and in- 
stanton (INST) character strings in JETSET 
common block 




LUEXEC 


Simulate the whole fragmentation and decay 
chain 




LULIST(l) 


List the (first) event record 




LUHEPC(l) 


Convert JETSET event record contents back 
to the HEPEVT common block 



Table 3 

Depending on the flag QICONT(18) = .TRUE. /.FALSE., HERWIG [8] and JET- 
SET [18] hadronization may be selected. The respective calls to HERWIG and 
(modified) JETSET routines used in QCDGEN are displayed. 
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QINWAR(l): 
QILPAR(8): 

QIUPAR(8): 
QIWARF(2): 



QINWAR(2): 
QILPAR(9): 

QIUPAR(9): 
QIWARF(3): 



QINWAR(3): 
QILPAR(6): 

QIUPAR(6): 
QIWARF(4): 



QIWARF(l) is reset to .TRUE.. 

= (default); counts the number of violations 

of Eq. (23). 

= ln(5 max (default); reset to the logarithm of 
the right hand side of Eq. (23) if this kinemat- 
ical constraint is violated. 
= <3 max ; reset to the right hand side of Eq. (23) 
if this kinematical constraint is violated. 
= .FALSE, (default); if Q'£ iD violates the kine- 
matical constraint 



Q 



r2 



> 



*^Bj min 



^max *^Bj min 



-wl 



(24) 



QIWARF(2) is reset to .TRUE.. 

= (default); counts the number of violations 

of Eq. (24). 

= ln<5 min (default); reset to the logarithm of 
the right hand side of Eq. (24) if this kinemat- 
ical constraint is violated. 
= <3 min ; reset to the right hand side of Eq. (24) 
if this kinematical constraint is violated. 
= .FALSE, (default); if x' mSuX violates the kine- 
matical constraint 



Wf ■ 

I < 1 _ £ mln 

•^max — 1 c 
^max*J 



(25) 



QIWARF(3) is reset to .TRUE.. 

= (default); counts the number of violations 

of Eq. (25). 

= lnx max (default); reset to the logarithm of 
the right hand side of Eq. (25) if this kinemat- 
ical constraint is violated. 
= x max ; reset to the right hand side of Eq. (25) 
if this kinematical constraint is violated. 
= .FALSE, (default); if x' min violates the kine- 
matical constraint 



•^min — 



(26) 



QIWARF(4) is reset to .TRUE.. 
QINWAR(4): = (default); counts the number of violations 
of Eq. (26). 

QILPAR(7): = In x' min (default); reset to the logarithm of the 
right hand side of Eq. (26) if this kinematical 
constraint is violated. 
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QIUPAR(7): 
QIWARF(ll): 

QINWAR(ll): 

QIC0NT(6): 

QIQMASS(1,I): 

QIQMASS(2,I): 
QIRCAL: 



= x' min ; reset to the right hand side of Eq. (26) 
if this kinematical constraint is violated. 
= .FALSE, (default); reset to .TRUE, if current 
quark mass requirement (c.f. QIC0NT(6)) does 
not fit to the selected order of variable genera- 
tion (c.f. QICONT(5)). 

= (default); counts the number of recurring 
QIWARF(ll) = .TRUE, settings. 
= .TRUE, (default); reset to .FALSE, if the flag 
QIWARF(ll) = .TRUE.. 

J2i=i f^qi-, cumulative sum of HERWIG [8] quark 
masses with 1 < I < 6 referring to the HER- 
WIG identity code. 

m qi ; HERWIG quark mass with 1 < I < 6 re- 
ferring to the HERWIG identity code. 
= .TRUE, (default); reset to .FALSE, on exit. 



Remarks: This routine is only performed if the flag QIRCAL is set to 
TRUE. (QIINIT, QISETD, QISETF, QISETI). 

SUBROUTINE QICCON 

Purpose: Assignment of the colour and flavour connections for the partons. 



Non-local variables (re)set: 

QIPLIS (JLP, ILP): IHEP pointer referring to the outgoing par- 
ton JLP in the string ILP. 
JMOHEP: array of HERWIG [8] "mother" pointers. 

JDAHEP: array of HERWIG "daughter" pointers. 



Remarks: For an explicit example, see the description of routine QISTID. 

Procedure: 

• The IHEP pointer of the current quark k (IHEP=10), as assigned in QI- 
HGEN, is added to the array QIPLIS of IHEP pointers at the position 
(JLP,ILP) of the virtual quark q'. 

• The IHEP pointer of the incoming gluon (IHEP=6), gned in QI- 
HGEN, is added to the array QIPLIS of IHEP pointers at the position 
(JLP,ILP) of the incoming gluon. 

• The colour flow is constructed by connecting the colour lines of neighbouring 
partons within each of the strings constructed in the subroutines QIGLST 
and QIGPAR. 

• Finally, the flavour flow is constructed by connecting the flavour lines of the 
quark at the beginning of a string with the flavour line of the anti-quark at 
the end of a string. 
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SUBROUTINE QIGETD (NUM, VALUE, OK) 

Purpose: Get value of double precision parameter QIUPAR(NUM) for given 
NUM. 

Arguments: 

NUM: integer, pointing to the parameter QIUPAR(NUM) to 
be read. 

VALUE: value of double precision parameter QIUPAR(NUM). 
OK: logical return flag. 

Non-local variables (re)set: 

QINWAR(IO): = QINWAR(10)+1; counts attempts to read non- 
existing parameter QIUPAR(.). 
QIWARF(IO): = .TRUE, if parameter index NUM out of range. 
Remarks: Service routine for reading out the double precision input param- 
eters QIUPAR(NUM), NUM=1,2,...,18. On return, the logical 
flag OK is .TRUE, if NUM is within its allowed range. 

SUBROUTINE QIGETF (NUM, FLAG, OK) 

Purpose: Get value of logical flag QICONT(NUM) for given NUM. 
Arguments: 

NUM: integer, pointing to the parameter QICONT(NUM) to 
be read. 

FLAG: boolean value of QICONT(NUM). 
OK: logical return flag. 

Non-local variables (re)set: 

QINWAR(6): = QINWAR(6)+1; counts attempts to read non- 
existing parameter QICONT(.). 
QIWARF(6): = .TRUE, if parameter index NUM out of range. 
QIRCAL: = .TRUE, if QIWARF(6) = .FALSE, (default). 

Remarks: Service routine for reading out the boolean input parameters QI- 
CONT(NUM), NUM=1,2,. . . ,21. On return, the logical flag OK 
is .TRUE, if NUM is within its allowed range. 

SUBROUTINE QIGETI (NUM, VALUE, OK) 

Purpose: Get value of integer parameter QIPARI(NUM) for given NUM. 



Arguments: 

NUM: integer, pointing to the parameter QIPARI(NUM) to 
be read. 

VALUE: value of QIPARI(NUM). 
OK: logical return flag. 
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Non-local variables (re) set: 

QINWAR(8): = QINWAR(8)+1; counts attempts to read non- 
existing parameter QIPARI(.). 

QIWARF(8): = .TRUE, if parameter index NUM out of range. 

QIRCAL: = .TRUE, if QIWARF(8) = .FALSE, (default). 

Remarks: Service routine for reading out the integer input parameters QI- 
PARI(NUM), NUM=1,2,. . . ,10. On return, the logical flag OK is 
.TRUE, if NUM is within its allowed range. 

SUBROUTINE QIGLST 

Purpose: Generation of n/ [q . . . g ] - "strings" of partons, with a total num- 
ber rig+1 of gluons inserted in between, in accord with the require- 
ment of "flavour democracy" among the rif gg-pairs participating 
in the instanton-induced subprocess. 

Non-local variables (re) set: 

QIPTYP(JLP,ILP): particle data group [16] identity code (ID- 

PDG) of the parton JLP in the string ILP. 

QINLIS (ILP): number of partons in the string ILP. 

QIPINF(JLP,ILP): = .FALSE.; logical flag indicating whether 

parton JLP in the string ILP 
is incoming. 

Procedure: 

Each of the rif generated strings begins with a quark, followed by a random 
number < n g + 1 of gluons, and ends with an anti-quark of randomly chosen 
flavour. Due to the required "flavour democracy", gg-pairs of all nf flavours 
must occur precisely once, e.g. for nf = 3 and n g = 3, 

[dgu] [uggs] [sgd]. (27) 

In practice, the strings are generated and administrated as follows: 

• First, a 2-dimensional integer array QIPTYP is defined, whose entries QIP- 
TYP (JLP, ILP) will contain the particle data group identity code (IDPDG) 
of the parton JLP in the ILP-th [q . . . q] - string. All entries are initially set 
to zero. 

• The quarks, labelled by their IDPDG-codes 1 ..rif, are then inserted con- 
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secutively at the top of the first nf columns ("strings"), e. g. for nf — 3, 





TT V> 

12 3 4 


30 




1 


1 2 3 0. 







2 


0. 





JLP 










30 


0. 





QINLIS(ILP) 


1110. 






• Next, the anti-quarks with IDPDG-codes —1 .. —rif, are inserted randomly 
in the second row of the first rif columns, e.g. 





1 


2 


ILP 

3 


4 


30 




1 


1 


2 


3 










2 


-2 


-3 


-1 








JLP 


3 


















30 

















QINLIS(ILP) 


2 


2 


2 









• Finally, the n g + 1 gluons of the instanton subprocess with IDPDG-code 21 
are distributed randomly over the first nf columns in between the quarks 
and anti-quarks: For each of the gluons, one of the first nf columns is picked 
at random, the anti-quark in this column is moved one row further down 
and the gluon is inserted at the anti-quark's previous position. A possible 
outcome of this procedure for the example above with n g = 3 is 





1 


2 


ILP 

3 


4 


30 




1 


1 


2 


3 










2 


21 


21 


21 










3 


-2 


21 


-1 







JLP 


4 





-3 












5 


















30 

















QINLIS(ILP) 


3 


4 


3 









SUBROUTINE QIGMUL (WGTFCT) 
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Purpose: Generation of the number of gluons n g emitted by the instanton 
subprocess in accord with kinematical constraints. 

Arguments: 

WGTFCT: = 1.0, returned, if the event is to be accepted. 

= 0.0, returned, if the final state quark masses 
are too high; event will be killed (default). 

Non-local variables (re)set: 

QINGLU: n g + 1, number of gluons involved in the instan- 

ton subprocess. 

QIWARF(14): = .FALSE, (default); reset to .TRUE, if the cal- 
culation of the number of gluons requires too 
many iterations. 

QINWAR(14): = (default); counts the number of recurring 
QIWARF(14) = .TRUE, settings. 

Procedure: 

• The minimal CM energy available for gluon generation is calculated and 
WGTFCT is set accordingly. 

• If the emission of at least one gluon is energetically allowed, the aver- 
age gluon multiplicity (n g )^ is calculated and returned by the function 
GMULT. 

• The number n g of emitted gluons is then generated according to a Poisson 
distribution [9] with mean (n g )^\ 




• An upper limit on (n g )^ is enforced and the kinematical constraints are 
checked again. In case they are violated, a new iteration cycle is started 
corresponding to a different generated value of n g . 

SUBROUTINE QIGPAR 

Purpose: Identification of the incoming partons among the partons partic- 
ipating in the instanton-induced subprocess. 

Non-local variables (re)set: 
QIPINC(1,1) 
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row number (JLP) in the 2-dimensional ar- 
ray QIPTYP (QIGLST), where the anti- 
parton corresponding to the incoming vir- 
tual (anti-) quark q' is located. 



QIPINC(1,2): 



QIPINC(2,1): 



QIPINC(2,2): 



QIPINF(JLP,ILP): 



column number (ILP) in the 2- dimensional 
array QIPTYP (QIGLST), where the anti- 
parton corresponding to the incoming vir- 
tual (anti-) quark q' is located, 
row number (JLP) of the incoming on-shell 
gluon g in the 2-dimensional array QIPTYP 
(QIGLST). 

column number (ILP) of the incoming on- 
shell gluon g in the 2-dimensional array QIP- 
TYP (QIGLST). 

= . FALSE, (default); set to .TRUE, only if 
(JLP, ILP) = (QIPINC(I,1),QIPINC(I,2)), 
with 1=1,2, corresponding to the incoming 
partons. 



Procedure: 

One of the n g + 1 gluons in the generated rif [q . . . g]-strings (QIGPAR) is 
randomly chosen and marked as incoming. Furthermore, the anti-parton cor- 
responding to the virtual (anti-)quark q' is identified and marked as incoming. 

SUBROUTINE QIHGEN 

Purpose: Main (hard) instanton process generator. 



Non-local variables (re) set: 

QISGAM: 

IDN(l): 

IDN(2): 

EMSCA: 

QIZPAR: 

QIXONE: 

QIYONE: 

QIXBJG: 

QIYBJG: 

QIQ2GA: 

QIELEN: 

QIPREN: 

EVWGT: 

QIMEWT: 

QIPINC(1,3): 



eP CM energy squared S. 
HERWIG identity code (IDHW) of incoming e. 
HERWIG identity code (IDHW) of incoming g. 
factorization scale fj,f — QIUPAR(18) 
generated z. 



x' z 



Q' 2 /{Sx'z) 
generated xbj- 
generated t/bj- 

Bjorken variable Q 2 calculated via Q 2 = SxBjl/By 
energy of incoming e, e t . 
energy of incoming P, P t . 

weight Wgp (39) corresponding to the instanton- 
induced cross section (3). 

internal storage of weight corresponding to Eq. (3) 
for analysis 

particle data group identity code (IDPDG) of 
virtual quark q'. 
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QIPINC(2,3): particle data group identity code (IDPDG) of 

incoming gluon g. 
QIPINC(1,4): (=10) IHEP pointer of current quark k. 
QIPINC(2,4): (=6) IHEP pointer of incoming g. 
QISHEP: CM energy squared of instanton subprocess, W 2 = 

(q> + g f = s >. 

JMOHEP: array of HERWIG [8] "mother" pointers. 

JDAHEP: array of HERWIG "daughter" pointers. 



Procedure: 



For the generation of an event, this routine is called twice: in the first call, a 
Monte Carlo weight associated with the instanton-induced total cross section 
is calculated. In the second call, the generation of the event is completed. 

In the first call (GENEV = .FALSE.), the following steps are performed in 
QIHGEN: 

• Information on the incoming beams is accumulated, initial state radiation 
from the lepton is optionally accounted for (EXFRAC), various derived 
parameters are calculated and kinematical limits are checked (QICALC). 

• The HERWIG [8] identity code (IDHW) of the current quark k and of the 
virtual quark q' (c. f. Fig. 1) are generated (QIHPAR). 

• Next follows the generation of the various Bjorken variables Q' 2 , x' ', z, xbj, 
l/Bj (c. f. Fig. 1) of the instanton-induced DIS process, with weights corre- 
sponding to the associated factors appearing in the instanton-induced total 
cross section (3): 

o By a call to QIHINS, the Bjorken variables Q' 2 and x' of the instanton 
subprocess are generated. Moreover, the Monte Carlo weight WQi x i a from 
Eq. (43), corresponding to the integral 

O' 2 x' 

max max j f n 

I dQ ' I f-i <$(*'<n (29) 

Q' 2 . x' • 

^mm mm 

is calculated, where o$ n denotes the instanton-induced q'q total cross sec- 
tion [4] from Eq. (17). At this stage, also the gluon multiplicity depending 
on x' and Q' 2 is generated (QIGMUL). 
o The fractional momentum z of the incoming gluon is generated and the 
Monte Carlo weight W zg corresponding to the integral 

2max 

J dzf g (z,(if), (30) 

( Q' 2 ^BjminA 
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is calculated, where f g is the gluon distribution in the proton. By default, 
f g is taken from Ref. [25] (set 1.1) and the factorization scale ///is identi- 
fied with Q' min . For the default setting of the control flags in QIINIT, the 
value of z as well as the Monte Carlo weight W zg are generated as follows: 
Random points y z , with 

y z = lnz <^> z = expy z (31) 

are generated uniformly in the interval between In zl and Inzjj, where zl 
and zjj are the lower and upper integration limits in Eq. (30), respectively. 
Hence, the Monte Carlo weight associated with Eq. (30) reads 

W zg = (\nzu -\nz L )z f g (z,/j, f ). (32) 

o The Bjorken variables xb] and yej are generated analogously to the z 
variable and the weight W Xw VBj , 

W XBi m = In i n 9(Sx Bm - Q 2 min ) , (33) 
XB]L Vb s l 

corresponding to the integral 

I dx^ I ^ Q{ _ QU (34) 

J xb\ J Ub\ 

is calculated. The integration limits in Eq. (34) are given in terms of the 
physical x B j, yej-cuts, with default values as displayed in Table 4, 

_ i m l 1 

X ~B] U — X Z q QI2~ j 

2/Bj max — 5^77 

•^Bj L *^Bj mini 

2/Bjt/ = 2/Bjmax; (35) 

( Q' 2 m 2 k 1 \ 

Finally, given the generated Bjorken variables, the weights associated with 
the remaining factors in the instanton-induced cross section (3) are calcu- 
lated: 



o The weight arising from the splitting of the photon into a gg-pair is ob- 
tained from the routine QISPLT, 

wv = ~t~ *#V» Q' 2 > *, Q 2 = Sx Bm )- (36) 

X z H 

The factor P^P from Eq. (4) accounts for the flux of virtual (anti-) quarks q' 
in the instanton background, entering the instanton-induced </g-subprocess 
from the photon side [4,14] (c.f. Fig. 1). 
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o The familiar weight from the flux of virtual photons is calculated, 

^-s^i 1 -**^)- (37) 

o The left-over weight factor required for a properly normalized cross section 
(in units of nb) is obtained from the routine QIPVWT, 

W n = 2 7T a 2 x Bj x' 4 • ( 38 ) 

q'=d,u,s,...; d,u,s,... 

Here a is the electro-magnetic fine-structure constant and e q < denotes the 
fractional electro-magnetic charge of the virtual (anti-)quark q'. The sum 
extends over quarks and anti-quarks of the considered rif light flavours. 

• The total event weight (EVWGT) required by the HERWIG package, is the 
product of all the above weights, 

= W Q , x , a W zg W XBi yBi W q , W, W n . (39) 

This ends the first call of QIHGEN. 

For standard settings of the control flags in the HERWIG [8] routine HWI- 
GIN (NOWGT=.TRUE.), the (standard) rejection method is applied pro- 
ducing an unweighted event distribution according to the instanton-induced, 
normalized differential cross section (9). 



If the event is accepted (GENEV = .TRUE.), QIHGEN is called again. The 
following steps are then performed: 

• The lepton beam e is stored in the event record. 

• The 4-momentum g of the incoming gluon is calculated from its gener- 
ated fractional momentum z and the 4-momentum of the incoming proton 
(QIKPAR). The result is stored in the event record. 

• The 4-momentum q of the virtual photon is generated (QIKGAM). 

• The 4-momentum e' of the outgoing lepton is calculated from e, q and stored. 

• The 4-momenta of the virtual quark q 1 and the outgoing current quark k 
are generated (QIKGSP) and copied into the event record. 

• The required data for the generation routines of the instanton subprocess 

q'g 4 X (40) 

are set up and its (total) momentum vector q' + g is calculated. 

• The colour connections of the hard subprocess are initiated. 

• The instanton-induced partonic final state X in Eq. (40) is generated (QIS- 
TID) with all colour connections in the g'g-CM system and, thereafter, 
transformed into the laboratory frame. 
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SUBROUTINE QIHINS 

Purpose: Generation of the Bjorken variables Q' 2 and x 1 ', associated with 
the instanton subprocess, and calculation of the Monte Carlo 
weight associated with the integral (29). 

Non-local variables (re) set: 

QIL2IN: Q' 2 according to Eq. (41). 

QIX1HT: x' according to Eq. (41). 

QIWGT: weight (43) associated with the integral (29). 

Procedure: 

The theoretical distributions in Q' 2 and x' are very strongly varying functions. 
Therefore, the improved generation strategy outlined in the description of the 
routine QIRDIS is used for default settings of the control flags in QIINIT. 

The variables Q' 2 and x' are generated by calls to the subroutine QIRDIS, 
which provides values of 

Q' 2 = {Z Q a)- l ' N Q'* and x' = (Z^)' 1 ^' , (41) 



with Zqi2 and Z x > uniformly distributed in the intervals 

< Zq>* < \ N „„ and - \ Ni < Z x , < 7-^r7> (42) 



respectively. The optimal powers Nqn (L2NV) and N x > (XPNV) are empiri- 
cally determined and set in QIHINS. 

In analogy to Eq. (67), the total weight WQ> x > a , associated with the integral 
(29), is given as a product of three weight factors 

W (yx , IT = W Q *W x ,W a , (43) 



where 



W r ,= 



1 



N Q ,2 
1 



(Q 



12 \—N r 



- (Q 



12 y N c 

max/ 



(Q>2yN Ql2 . 



r > )-N x , _ ( r > )-N x 



1 



(X>) 



W^^Q' 2 a^{x\Q' 2 ). 



x 



(44) 
(45) 
(46) 



The function Q2SIG is invoked to calculate the weight factor (46) from the 
instanton-induced q'g total cross section. 
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SUBROUTINE QIHPAR 



Purpose: Generation of the HERWIG [8] identity codes (IDHW) of the cur- 
rent quark k as well as of the virtual quark q'. 

Non-local variables (re)set: 

QINFAM: = QIPARI(4); number of (light) flavours n f . 

QIPINT(l): IDHW of current quark k. 

QIPINT(2): IDHW of virtual quark q 1 . 

QICQM2: m\\ HERWIG mass squared of current quark. 

Procedure: 

By means of the rejection method, the identity code IDHW of the current 
quark or anti-quark k is chosen randomly among 2n/ light quarks and anti- 
quarks, {d, u, s, . . . ; d,u,s,...}, according to the probability distribution 

e 2 

— ^ -; i = d,u,...;d,u,..., (47) 

q=d,u,...; d,u,... 



reflecting the flavour structure of the standard electromagnetic 7% (^-coupling. 
The identity code of the virtual quark q' is then given in terms of IDHW(/c). 
In Eq. (47), e; is the electric charge of the quark flavour i. 

SUBROUTINE QIINIH 

Purpose: Initialization of "particle" 'INST' in HERWIG [8] event record. 

Procedure: Association of the name 'INST' with the HERWIG identity code 
(IDHW) '206' of the instanton "particle" in the HERWIG event 
record. 

SUBROUTINE QIINIT 

Purpose: Initialization of input parameters as displayed in Tables 4 and 5. 

Procedure: 

The routine is designed to make the changing of default values as easy as pos- 
sible. The parameters available to the user are listed in Tables 4 and 5 (for the 
definition of kinematic variables see Fig. 1). Furthermore, a name for the "in- 
stanton particle" is set by QIPARC(1)='INST'. After setting the parameters, 
QIINIT calls the subroutine QIINIH to initialize the new instanton particle 
in the HERWIG event record. Finally, it sets the flag QIRCAL = .TRUE. (<£> 
calculate derived quantities and check limits). 
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Name 


Default value 


Description 


QIUPAR(l) 


0.75 


Effective gluon mass 


QIUPAR(2) 


1.0 • 10" 3 


^Bjmin: minimum allowed XBj 


QIUPAR(3) 


1.0 


z max : maximum of momentum 
fraction z carried by the gluon 


QIUPAR(4) 


1.0 


2/Bjmax: maximum allowed i/Bj 




n 1 
0.1 


yBjmin : minimum allowed i/Bj 


QIUPAR(6) 


0.9 


x' mSiX : maximum allowed x' 


QIUPAR(7) 


0.35 


£min : minimum allowed x' 


QIUPAR(8) 


(v/QIUPAR(9) + 30) 2 


Qmax : maximum allowed Q' 2 


QIUPAR(9) 


( Q QIUPAR(13) \ 2 
[° LAMNF (3, 0.15267) J 


Qmm'- rninimum allowed Q' 2 


QIUPAR(IO) 


QIUPAR(7) 


Cut on x' below which the cross 
section is assumed constant 


QIUPAR(ll) 


-1.0- 10~ 10 


Minimum allowed weight for exit 

Ul Illcllll WclgilL gcllcl dllUIl 


QIUPAR(12) 


0.0 


Wi m - m : minimum allowed instan- 
ton CM energy 


QIUPAR(13) 


LAMNF(QIPARI(4), 0.219) 


A^i : from input value (3 loop) 
A§ = 0.219±8;gg GeV [16] 


/ \ T T T|") A T"> / -j r\ 

(^lUrAK(lo) 


n n 

0.0 


Minimum energy after mass sub- 
traction from instanton 


QIUPAR(16) 


0.15 


(i r /Q' : renorm. scale in Q' units 


QIUPAR(17) 


QIUPAR(9) 


Q 2 nin : minimum allowed Q 2 


QIUPAR(18) 


V /QIUPAR(9) 


fij: factorization scale 


QIPARI(l) 


206 


HERWIG [8] identity code 
(IDHW) of instanton "particle" 


QIPARI(2) 


100 


Maximum number of MAMBO 
iterations per weight 


QIPARI(3) 


300 


Maximum number of phase space 
iterations per MAMBO weight 


QIPARI(4) 


3 


nf. number of (light) flavours 


QIPARI(6) 


10 


Maximum allowed average gluon 
multiplicity 


QIPARI(7) 


20 


Maximum number of iterations 
in the cross section weight gen- 
eration step 


QIPARI(9) 


40 


Maximum number of QIGMUL 
iterations 


QIPARI(IO) 


3 


Number of loops in a^g evalua- 
tion 


QINWAR(I) 





Counts recurrence of warning I, 
I = 1..14 (c.f. Table 5) 



Table 4 

Floating point and integer parameters set in QIINIT. All energy/mass dimensions 
are in GeV 
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Name 


Default 


Description 


QTCONTfFl 


T 


TTcp ti on-tri^n a 1 ViriPTnatip Txrpichf" tot omrpn Tihfitip git*sipp 

\J DC 11U11 LllVldl IVIIICIIIOjLIC WCl^jlllj 1U1 glVCll JJlldDC O^JCICC 






distribution 


QIC0NT(2) 


T 


Use standard routine (leading-order) to generate 






weight for given phase space distribution 


QIC0NT(3) 


T 


Use weight due to instanton cross section 


QIC0NT(4) 


T 


Disregard instanton minimum mass requirement 


QIC0NT(5) 


T 


Generate Q' 2 before x' 


QIC0NT(6) 


T 


Enforce mass of current quark in kinematics 


QIC0NT(9) 


T 


Enforce maximum allowed number of gluons 


QICONT(ll) 


T 


Check that energy suffices for gluon generation 


QIC0NT(13) 


T 


Kill events with mass too high 


QIC0NT(14) 


T 


Use z generation as dz/z 


QIC0NT(15) 


T 


Generate x' with efficiency parametrization 


QlCONT(16) 


I 


Generate Q with efficiency parametrization 


QIC0NT(18) 


T 


Use HERWIG rather than JETSET hadromzation 


QlCONT(20) 


I 


Azimutnal angle of e generated randomly 


QIC0NT(21) 


T 


run nadromzation on/ on 


QIWARF(l) 


F 


Q 12 upper generation limit reset 


0IWARFf2) 


F 


Q lower generation limit reset 


0IWARFf3) 


F 


x f upper generation limit reset 


v^l V V r\ 111 I J 




Ju lUWci ^ciiciaLiun 111I11L ICbCt 


v^l V V Jr\ 111 IO 1 




iT.L LCllip L LU llclDCt IIUII-CXIdLcIIL lU^lL-dl p<al cllllt; Lcl 








^1 X\ t\ 1 VI ^0 J 


r 


r\.i Leinp l lo i eau noii-exis LeiiL logical pai ameiei 






fOTdFTF"! 


Vc^l V \ -Til V 1 ^ ( J 


TP 


r\i Lcinp l lo yL c jbci iioii-exibiciiL iiiieg,ei pai ameiei 








OTWAT^ThYR^ 
^1 VVrvIxr ^0 j 


r 


rA.L Leinp l lo i cdu non-exiDLeiiL niLc^d pai amcL(ji 






(T)TnF,TT^ 




TP 

r 


A TT /"V TTT T~\ 1 T~ y^l / Y*{~\ \ C /Z1T TT /~T T^l AVI C""T /~\TT T" /H T 1 £~l T\ T* O 1 C 1 i^lT^l 7 "\ o 

ritxempi lo ^lejsei non-exisieiit oouoie piecihion pa- 






rameter fOISETD'l 


QIWARF(IO) 


F 


Attempt to read non-existent double precision param- 






eter (QIGETD) 


QIWARF(ll) 


F 


Current quark mass requirement does not fit to the se- 






lected order of variable generation (c. f. QIC0NT(5)) 


QIWARF(12) 


F 


Negative k\ for current quark, reset to 0.0 


QIWARF(14) 


F 


Too many iterations in QIGMUL 


QIRCAL 


T 


Check kinematical limits 



Table 5 

Logical flags set in QIINIT (T = .TRUE., F = .FALSE.). 
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SUBROUTINE QIKGAM 



Purpose: Generation of the 4-momentum q of the virtual photon, with the 
constraint —q 2 = SxbjVb] = Q 2 > 0. 



Non-local variables (re)set: 

QIPGAM (I): I-th component of the 4-momentum of the virtual 
photon, q = (q x ,q y ,q z ,q t ), where I={1,2,3} O 
{x, y, z) and 1=4 t. 

Procedure: 



With the convention e = (0, 0, e z , et), e z > 0, for the incoming lepton momen- 
tum in the laboratory frame and assuming e 2 = 0, the 4-momentum of the 
virtual photon is expressed as 



q x = ~ yJSxBjVBjil — 2/Bj) cos <f> q , (48) 
% = ~\] Sx B] y B - s {l - y B] ) sin0 g , (49) 

q z = yB- s e t + 3 J , (50) 
4e t 

* = (51) 

according to a Sudakov decomposition. The azimuthal angle (f> q , —n < 4> q < 
+7T, is randomly generated for the default setting of the control flag QI- 
CONT(20) = .TRUE.. Otherwise, <p q is fixed at zero, i.e. e' (the 4-momentum 
of the scattered e* 1 ) and q are lying in the x — z plane. 

SUBROUTINE QIKGSP 

Purpose: Generation of the 4-momenta q' and k of the virtual quark and 
the outgoing, on-shell current quark, respectively (c. f. Fig. 1). 
Implementation of the constraints — q' 2 = Q' 2 > 0, where Q' 2 is 
the virtuality generated in QIHINS, and k 2 = m\. 



Non-local variables (re)set: 

QIQGAM (I, 2): I-th component of the 4-momentum of the vir- 
tual quark, q' = (q' x , q' y , q' z , q' t ), where I={1, 2, 3} 
<^> {x, y, z] and 1=4 t. 

QIQGAM (1, 1): I-th component of the 4-momentum of the cur- 
rent quark, k = (k x , k y , k z , k t ), where I={1, 2, 3} 
<^ {x, y, z} and 1=4 <^ t. 

QIWARF(12): = .FALSE, (default); reset to .TRUE, if the 
discriminant D 2 < in Eq. (57). 
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QINWAR(12): = (default); counts the number of recurring 
QIWARF(12) = .TRUE, settings. 

Procedure: 

By convention, the laboratory frame corresponds to the incoming proton and 
the incoming lepton travelling in — z and +z direction, respectively. Using a 
Sudakov decomposition to incorporate the required contraints on q' 2 and k 2 , 

- q ' 2 = Q' 2 > and k 2 = m 2 k , (52) 



the 4-momentum q' of the virtual quark is expressed as 



Q' 2 i 

q' x = D cos</y - — JSxBjVBjil ~VBj) cos g , 

ox z y-Qj 

Q' 2 i 

q' y = D sin</y - — \] Sx Bj y Bj (l - y B j) sin0„ 

OX Z ?/Bj 



(53) 
(54) 



q' z = 2D (cos^/cos^ + sin^/sin^) J Sx Bi y B j(l - y Bj )- — (55) 



12 



Q 

bx'z 



Q 



12 



Sx'zy B] 



(X B] - x'z + X B j(l - 2/Bj)) - %Bj 



mi 



Sy 



q' t — —2D (cos (j) q > cos 4> q + sin <f> q , sin 9 ) ^Sx Bi y Bi (l - y B] ) 



Sy B j 



(56) 



Q' 2 n ( Q' 2 . , „ u \ 

+ Pt [s*^ (xBj - ^ + ^(i - y Bj )) - x Bi - —j , 



with 



D= (57) 

„ Q' 2 / q /2 a n/2 / g /2 \ g /2 2 

^ Bj?/Bj ^ Bj V^'^bj ) Q \Sx>zy Bi V 5^ym mfc - 

The azimuthal angle </y, — n < <fi q r < +7r, is generated randomly. The 4- 
momentum of the current quark is then calculated via 

k — q — q', (58) 



where the photon 4-momentum q has been generated in QIKGAM. 
SUBROUTINE QIKPAR 

Purpose: Calculation of the 4-momentum g of the incoming gluon from its 
generated fractional momentum z and the 4-momentum P of the 
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incoming proton. 

Non-local variables (re) set: 

QIPPAR (I): I-th component of the 4-momentum of the incom- 
ing gluon, g = (g x , g y , g z , g t ), where I={1, 2, 3} 
{x, y, z} and 1=4 <^> t. 

QIPPAR (5): m g ; HERWIG [8] mass of the incoming gluon. 

Procedure: 

The momentum vector g of the incoming, on-shell gluon (with effective (HER- 
WIG) mass ^fg 2 = m g ) is calculated in the laboratory frame according to 



9 = {9x, 9 y , 9z, 9t) = (59) 
-[0,0,z(P t + \Pz\)- 7WT1-B n ,z( p t+ \Pz\)+ " ?/ 1 



z(p t +\p z \y v 11 1 zu Z (p t +\p z \)j> 

where z denotes the generated fractional momentum of the incoming gluon 
and P = (0, 0, P z , P t ) is the 4-momentum of the incoming proton. If P z < 0, 
the sign of g z is reversed, g z — > —g z . 

SUBROUTINE QIPLST 

Purpose: Assignment of masses for the 2 rif — 1 + n g outgoing partons from 
the instanton-induced subprocess. 

Non-local variables (re)set: 

QIPHEP (5, J): mass of the J-th outgoing parton from the 
instanton-induced subprocess. 

QINHEP: 2nj — 1 + n g ; number of outgoing partons 

from instanton subprocess. 

Procedure: 

To each outgoing parton an appropriate mass (default HERWIG [8] mass for 
quarks/ ant i- quarks, default gluon mass from QIINIT) is assigned. 

SUBROUTINE QIPSGN 

Purpose: Generation of the 4-momenta for the outgoing partons from the 
instanton-induced subprocess. 

Non-local variables (re) set: 

QIPHEP (I, J): I-th component of the 4-momentum corre- 

sponding to the J-th outgoing parton from 
the instanton subprocess. Here, I = {1, 2, 3} 
^ {x, y, z} and I = 4 t. 
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Procedure: 



• By means of the MAMBO [26] algorithm^], the momenta pj of the n = 
2nf — 1 + n g outgoing partons in the q'g CMS are generated uniformly in 
phase space, 

I n {d% { P ) - mj) } 6® £ p k )S[ Wj - £ p{ • (60) 
i=i \fc=i / \ fe=i / 

The masses were assigned in the subroutine QIPLST. 

• Next, for the default setting of the control flags in QIINIT (QICONT(l) = 
.TRUE, in Table 5), different energy weights for gluons and quarks (QIP- 
SWT) are implemented by means of the rejection method. The resulting 
momenta pj are then uniformly distributed in energy-weighted phase space, 
Eq. (10), corresponding to the leading-order matrix elements [3]. 

SUBROUTINE QIPSTO 

Purpose: Store 4-momenta and particle identity codes of all outgoing par- 
tons in the PHEP common block of HERWIG [8] and set appro- 
priate IHEP pointers. 



Non-local variables (re)set: 

PHEP (I, IHEP) : I-th 4-momentum component of the J-th out- 
going parton from the instanton subprocess, 
with IHEP=10+J and 1 < J < 2n f - l+n g . 

IDHEP(IHEP): particle data group [16] identity code (ID- 

PDG) of the J-th outgoing parton. 

IDHW (IHEP) : HERWIG identity code of J-th outgoing par- 
ton. 

QIPLIS (JLP, ILP): IHEP pointer referring to the outgoing par- 
ton JLP in the string ILP. 
Remarks: See also the descriptions of the routines QISTID and QIGLST. 

SUBROUTINE QIPSWT (PSWGT) 

Purpose: Calculation of the relative energy-weight factor (c. f. Eq. (10)) for 
the generated momenta of the outgoing partons, as correspond- 
ing to the modulus squared of the leading-order matrix elements. 



Argument: PSWGT: relative energy weight (63). 
Procedure: 

5 The MAMBO algorithm implemented in the Monte Carlo generator for baryon 
number violating interactions HERBVI [23] is actually used and provided in the 
source file 'saddle. F' of the QCDINS distribution. 
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An energy- weight factor w corresponding to the modulus squared of the leading- 
order matrix elements [3] is calculated as follows^: Each outgoing quark with 
four- momentum p q is weighted by its energy p*, each outgoing gluon with 
four- momentum p k by its energy squared p k 2 , such that 



™= II p\ UpI 2 - (61) 

q=l k=l 



It is easy to show that the corresponding maximum weight is given by 

Wr 



w = 2 2 ™ 9 

"'max " 



2 (n g +nf) — 1 



2 {tig+rif)— 1 

(62) 



The subroutine QIPSWT returns the relative weight 

w ve \ = w/w max (63) 

in the variable PSWGT. 
SUBROUTINE QIPVWT (WT) 

Purpose: Calculation of the remaining weight factors required for a prop- 
erly normalized cross section (3). 

Argument: WT: returned weight factor (64). 
Procedure: 



The weight factor 



W n = 2 7ra 2 x Bj x' ]T ej (64) 



q 

q'=d,u,s,...; d,IT,s,... 



is calculated. Here, a is the electro-magnetic fine-structure constant and e q i de- 
notes the fractional electro-magnetic charge of the virtual (anti-)quark q' . The 
sum extends over quarks and anti-quarks of the considered n f light flavours. 

SUBROUTINE QIRDIS (LU, LL, X, WGT, N) 

Purpose: Strong increase of efficiency through generation of X = Q' 2 , x' as 
dX/X N+1 . 



Arguments: 

6 We thank M. Seymour for pointing out an error in the previous version of this 
subroutine and for providing us with a corrected version [27] of it. 
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LU: logarithm of upper integration limit Xu in Eq. (65). 
LL: logarithm of lower integration limit Xl in Eq. (65). 
N: power in Eq. (66). 

X: X value associated with the uniformly distributed Z ac- 
cording to Eq. (66). 
WGT: Monte Carlo weight (67) with / = 1. 
Remarks: The instanton-induced cross section (3) involves integrations 

HXu,X L )= J ^f(X), (65) 

x L 

where, typically, the function f(X) grows steeply towards small X. 
This makes the standard Monte Carlo evaluation (e.g. Eqs. (31), 
(32)) of the integral (65) very inefficient. 

Procedure: 

The efficiency is strongly increased by means of the present routine. First, the 
integration variable is changed to 

Z = X~ N ; dZ = -NX- (N+1) dX, (66) 



with the power N empirically optimized and typically ranging between 2.5 and 
5.0. Uniformly distributed random points Z = X~ N are then generated in the 
interval X^ N < Z < X£ N . With this procedure, the weight corresponding to 
the integral (65) takes the form 

W = ±(X E N - X V N ) /(X= /~ 1/jV) . (67) 
SUBROUTINE QISETD (NUM, VALUE, OK) 

Purpose: Set double precision parameter QIUPAR(NUM) to a desired value. 
Arguments: 

NUM: integer, pointing to the parameter QIUPAR(NUM) to 
be changed. 

VALUE: desired double precision value of QIUPAR(NUM). 
OK: logical return flag. 

Non-local variables (re)set: 

QINWAR(9): = QINWAR(9)+1; counts attempts to (re)set 

non-existing parameter QIUPAR(.). 
QIWARF(9): = .TRUE, if parameter index NUM out of range. 
QIRCAL: = .TRUE, if QIWARF(9) = .FALSE, (default). 
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Remarks: Service routine for resetting the default double precision input 
parameters QIUPAR(NUM), NUM=1,2,. . . ,18, set in QIINIT 
as in Table 4, to the desired value VALUE. On return, the logical 
flag OK is .TRUE, if NUM is within its allowed range. 

SUBROUTINE QISETF (NUM, FLAG, OK) 

Purpose: Set logical flag QICONT(NUM) to a desired value. 

Arguments: 

NUM: integer, pointing to the parameter QICONT(NUM) to 
be changed. 

FLAG: desired boolean value of QICONT(NUM). 
OK: logical return flag. 

Non-local variables (re) set: 

QINWAR(5): = QINWAR(5)+1; counts attempts to (re)set 

non-existing parameter QICONT(.). 
QIWARF(5): = .TRUE, if parameter index NUM out of range. 
QIRCAL: = .TRUE, if QIWARF(5) = FALSE, (default). 

Remarks: Service routine for resetting the default boolean input parameters 
QICONT(NUM), NUM=1,2,. . . ,21, set in QIINIT as in Table 5, 
to the desired value FLAG. On return, the logical flag OK is 
.TRUE, if NUM is within its allowed range. 

SUBROUTINE QISETI (NUM, VALUE, OK) 

Purpose: Set integer parameter QIPARI(NUM) to a desired value. 
Arguments: 

NUM: integer, pointing to the parameter QIPARI(NUM) to 
be changed. 

VALUE: desired integer value of QIPARI(NUM). 
OK: logical return flag. 

Non-local variables (re)set: 

QINWAR(7): = QINWAR(7)+1; counts attempts to (re)set 

non-existing parameter QIPARI(.). 
QIWARF(7): = .TRUE, if parameter index NUM out of range. 
QIRCAL: = .TRUE, if QIWARF(7) = FALSE, (default). 

Remarks: Service routine for resetting the default integer input parameters 
QIPARI(NUM), NUM=1,2,. . . ,10, set in QIINIT as in Table 4, 
to the desired value VALUE. On return, the logical flag OK is 
.TRUE, if NUM is within its allowed range. 
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SUBROUTINE QISPLT (WEIGHT) 

Purpose: Calculation of the weight associated with the flux of the virtual 
(anti-)quark q' from the 7*g'/c vertex (c. f. Fig. 1). 

Argument: WEIGHT: returned weight factor (68). 
Procedure: 

The weight 

W «' = j- z P ? (68) 

is calculated, where PjjP is the flux factor (4) from Refs. [4,14]. 
SUBROUTINE QISTAT 

Purpose: Print input parameter vectors QIUPAR, QIPARI and QICONT, 
as well as possible warnings (QIWARF). For default values, see 
Tables 4 and 5. 

SUBROUTINE QISTID 

Purpose: Generation of identity codes, 4-momenta and colour /flavour con- 
nections of the 2n,f — 1 + n g outgoing partons of the instanton- 

induced subprocess g + q' 4> (2 n/ — 1) ■ q + n g ■ g . 

Procedure: 

• Given are the identity codes of the incoming partons, the number n g + 1 
of gluons and the required "flavour democracy" of the possible partonic 
subprocesses like e. g. for rif = 3 and n g = 4, 

s + g J^d + d + u + u + s + 4- g. (69) 

The structure as well as the identity codes of the partonic final state are 
then set up by the following two steps: 

o First, rif [q. . .q] - "strings" of partons are generated (QIGLST). Each 
begins with a quark, followed by a random number < n g + 1 of gluons, 
and ends with an anti-quark of randomly chosen flavour. Due to flavour 
democracy, gg-pairs of all rif flavours must occur precisely once, e.g. for 
the process (69) above, 

[dggs] [ugd] [sggu\. (70) 

All partons are marked as outgoing at this stage. 
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o Next (QIGPAR), one of the n g + l gluons in the above strings is randomly 
marked as incoming. Furthermore, the anti-parton corresponding to the 
virtual (anti-)quark q' is identified and marked as incoming, e.g. (•) for 
the process (69) and string configuration (70) above, 

[dgg s] [ugd] [s g g u] . (71) 



• The assignment of the final state parton momenta is based on the crucial 
requirement of isotropy in the "instanton rest frame", q' + g = 0. It proceeds 
through the following steps: 

o The appropriate (HERWIG [8]) mass m 8 is assigned to each of the 2nf — 

1 + n g outgoing partons i (QIPLST). 
o The 4-momenta pi in the instanton rest frame are generated uniformly 

in energy-weighted phase-space (10), corresponding to the leading-order 

matrix element [3] with different energy weights for gluons and quarks 

(QIPSGN). 

o The four-momenta and particle identity codes of all outgoing partons are 
stored in the PHEP common block of HERWIG (QIPSTO). The corre- 
sponding IHEP pointers start at IHEP=11 for the first outgoing parton 
from the instanton subprocess (a d quark in our example (71) above) and 
are increased consecutively along the strings, e. g. for the process (69) and 
string configuration (71) above, the IHEP pointers would read 

• # 

[ d g g s ] [ugd] [ s g g u ] 

11 12 13 14 15 16 17 18 19 . 

• Finally (QICCON), the stored list of IHEP pointers is augmented by those 
of the current quark (IHEP=10) and of the incoming gluon (IHEP=6), as 
assigned in QIHGEN. For the example above, this means 

• # 

[ d g g s ] [ugd] [ s g g u ] 

11 12 13 10 14 15 16 17 18 6 19 . 

As a main issue, the colour /flavour connections are established in QIC- 
CON as follows: The colour flow is constructed by connecting the colour 
lines of neighbouring partons within each of the strings. The flavour flow is 
constructed by connecting the flavour lines of the quark at the beginning of 
a string with the flavour line of the anti-quark at the end of a string. For 
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IHEP 


Description 


JMOHEP(2,IHEP) 


JDAHEP(2,IHEP) 


11 


d out 


12 


10 


12 


g out 


13 


11 


13 


g out 


10 


12 


10 


~s out (current quark) 


11 


13 


14 


7/ out, 


15 


16 


15 


out 


16 


14 


16 


<i out 


14 


15 


17 


s out 


18 


19 


18 


g out 


6 


17 


6 


5 in 


19 


18 


19 


¥ out 


17 


6 



Table 6 

HERWIG [8] mother /daughter pointers specifying the colour /flavour flow for the 
example (69) in the description of subroutine QISTID. 

the example above, the corresponding HERWIG mother/daughter pointers 
specifying the colour /flavour flow are then as in Table 6. 

SUBROUTINE QIUSPS (WT) 

Purpose: Dummy routine for optional modification of relative phase space 
weight as calculated in QIPSWT. 

Arguments: WT: modified relative phase space weight. 

Remarks: Called from QIPSWT if QICONT(2) = .FALSE.. Calculation of 
WT has to be provided by the user. 

FUNCTION XI (XPR, XLMIN, XLMAX, XQ, XMU, NF, LOOPFL) 

Purpose: The saddle point value of the conformally invariant //-separation 
£ is calculated as function of x', An / OLy^{Q') , An / a^(/i r ) and rif. 

Arguments: 



XPR: x' 

XLMIN: £ min ; lower boundary of £ used for interpolation. 

XLMAX: £ max ; upper boundary of £ used for interpolation. 

XQ: 4n/a m (Q>) 

XMU: 4ir/a m (fi r ) 

NF: ri/; number of light flavours. 

LOOPFL: =1: 1-loop renormalization group (RG) invariance [4] 
along with 1-loop form of a s . 
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= 2: 2-loop RG invariance [4] along with 2- loop form 
of a m . 

= 3: (default) 2-loop RG invariance along with 3-loop 

Procedure: 



form of a^As- 



First, the saddle point equation (18) is solved analytically for x', and the 
x--values corresponding to a discrete set of £j values are calculated as 

x < = j& ~ 2 ) (72) 

1 (6 + 2)+45(6)(5(6)-2« i )/T;?' 



with £ min < ^ < £ max and Vi from Eqs. (19)-(22) inserted: 



v t = v l (^X(Q'),X( f i r )) = 2D(S) x 



(73) 



exp{|- 



+ 



1 Ai-l 



(5) 1 2 Ai/3 



At In f AlX(Mr) ^ - ll X(u ) + 



V 



2,0(5) 



/ 



In Eq. (73), the explicit dependence of Vi on Q'//i r (c. f. Eq. (19) has been 
eliminated in favour of an X(Q'), X(/i r ) dependence by means of the standard 
scale transformation of a^g. Here, we have introduced the shorthand 



X(M) = 



An 
«Ms( M )' 



(74) 



The desired continuous inversion = £*(£', . . .) is then achieved by means of 
numerical interpolation based on the above exact supporting points 

FUNCTION XQS (QLAM, LOOPFL, NF) 

Purpose: Calculation of Aix/a-y^ as function of M/A^\ where M is a 
generic mass scale. 



Arguments: 

QLAM: M/Ag| } 

LOOPFL: = 1: 1-loop form of a s . 

= 2: 2-loop form of ct^g. 

= 3: (default) 3-loop form of a^g- 
NF: rif] number of light flavours. 
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Procedure: 



The running coupling is calculated according to the explicit formula 
Eq. (9.5a) in Ref. [16], which is accurate to 3-loop. The quantity 47r/«j^g 
is returned. The integer variable LOOPFL = 1, 2, 3 controls the loop-order 
at which Eq. (9.5a) in Ref. [16] is truncated. 



4 Usage and availabiliy 



QCDINS 2.0 should be loaded together with HERWIG [8] version 5.9 and 
JETSET [18] version 7.4, that are part of the CERNLIB distribution. The 
program is a slave system, which the user must call from his own steering pro- 
gram. A very simple example is provided in Appendix A. By default, QCDINS 
is compiled in form of a library, 'libqcdins.a', that may be linked together with 
'libherwig59.a', 'libjetset74.a' and the library 'libpdflib.a' of parton distribu- 
tion functions to the steering program. An extensive demonstration program, 
including an interface to the HZTOOL [28] package, and a detailed installation 
instruction are included in the distribution. 

Information about the distribution, its update history, an interactive manual, 
the source code, pictures of typical events etc. can be accessed via the QCDINS 
WWW site, |http: / / www.desy.de/~t00fri / qcdins / qcdins.html . 
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Appendix A: Simple steering program 



The following code may serve as the simplest illustration of a steering pro- 
gram for the QCDINS/HERWIG package. Note that the two subroutines 
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HWANAL and HWAEND must exist. 



For an example of an interface with the JETSET hadronization, we refer to 
the extensive steering program 'qtesthz' included in the QCDINS distribution. 



PROGRAM QCDINS 
#include "herwig.inc" 
C Force inclusion of block data 

EXTERNAL HWUDAT 

C Initialize process number 

IPROC = 17600 
C Maximal number of events in this run 

MAXEV = 1000 
C Beam particles 

PBEAM1 = 27.5D0 

PBEAM2 = 820. 0D0 

PARTI = >E+' 

PART2 = >P> 
C Initialize common blocks 

CALL HWIGIN 

C User can reset parameters at this point, otherwise values 
C set in HWIGIN will be used. 
C No vertex information in event printout 
PRVTX=. FALSE. 

C Reset number of shots for initial max weight search 

IBSH = 5000 

LRSUD=0 

LWSUD=77 
C Seeds 

NRN(1)=106645412 

NRN(2)=135135176 
C Use laboratory frame 

USECMF = .FALSE. 
C Compute parameter-dependent constants 

CALL HWUINC 
C Number of HERWIG events to print out 

MAXPR = 1 

C Call hwusta to make any particle stable 
CALL HWUSTA ('PI0 ') 
CALL HWUSTA ('K_S0 ') 

C Initialize default QCDINS input parameters 

CALL QIINIT 
C Print input parameters 

CALL Q I STAT 
C Loop over events 

CALL QCL00P 

C User's terminal calculations 
CALL HWAEND 
STOP 
END 

SUBROUTINE HWANAL 

RETURN 

END 

SUBROUTINE HWAEND 

RETURN 

END 
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Appendix B: Test run output 



Below, we display the essential output from a test run of the very simple 
steering program from Appendix A. This test run simulates 1000 complete 
instanton-induced events in deep-inelastic e + P scattering (HERA), with E e + = 
27.5 GeV and E P = 820 GeV, in the laboratory frame. All the QCDINS pa- 
rameters correspond to the default values as set in the initialization routine 
QIINIT (c.f. Tables 4 and 5). 

For reasons of space, the display of the first event has been truncated after 
the hard subprocess level. 



INPUT CONDITIONS FOR THIS RUN 

BEAM 1 (E+ ) MOM. = 27.50 

BEAM 2 (P ) MOM. = 820.00 

PROCESS CODE (IPR0C) = 17600 



QCD Instanton Monte Carlo Information Version 2.0 



Parameter Value 



Default gluon mass (GeV) 

Minimum allowed value of x_bj 

Maximum allowed value of z (proton mom.frac.) 

Maximum allowed value of y_bj 

Minimum allowed value of y_bj 

Maximum allowed value of x prime 

Minimum allowed value of x prime 

Maximum allowed value of Q prime **2 (GeV**2) 

Minimum allowed value of Q prime **2 (GeV**2) 

Lower cut for ME calculation on x prime 

Lowest allowed weight efficiency cut 

Minimum instanton invariant mass 

Lambda-MS-bar(nf) [GeV] from PDG 1998 

Minimum total K.E. of outgoing partons (GeV) 

Renormalization point KAPPA = mu_r/Qprime : 

Minimum allowed value of Q **2 (GeV**2) 

Factorization scale (GeV) 



0. 


7500E+00 


0. 


. 1000E-02 


0. 


1000E+01 





. 1000E+01 





1000E+00 





9000E+00 


0. 


3500E+00 


0. 


1652E+04 


0. 


1134E+03 


0. 


3500E+00 




. 1000E-09 


0. 


0000E+00 


0. 


3459E+00 


0. 


0000E+00 


0. 


1500E+00 


0. 


1134E+03 


0. 


. 1065E+02 



ID code assumed for instanton 206 

Maximum number of MAMB0 iterations per PS wt . 100 

Maximum number of PS wt . rejections per event 300 

Number nf of (light) flavours 3 

Maximum average gluon multiplicity in distbn. 10 

Maximum number of iterations for ME generation 20 

Maximum number of QIGMUL iterations 40 

Number of loops in RG-invariance/alpha_s 3 
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Control flag option 


Setting 


Reweight phase space configurations 


True 


Use default phase space reweighting 


True 


Use Matrix element weight 


True 


Disregard instanton minimum mass requirement 


True 


Generate Q prime before x prime 


True 


Enforce mass of current quark in kinematics 


True 


Enforce limit on maximum number of gluons 


True 


Ensure mass less than subprocess energy 


True 


Kill events with insufficient instanton mass 


True 


Use z generation as dz/z 


True 


Use x prime **-n generation for efficiency 


True 


Use Q prime **-n generation for efficiency 


True 


Use HERWIG rather than JETSET hadronization 


True 


Use random azimuth angle for scattered electron 


True 


Use full hadronization 


True 



INITIAL SEARCH FOR MAX WEIGHT 

PROCESS CODE IPROC = 17600 

RANDOM NO. SEED 1 = 1246579 

SEED 2 = 8447766 

NUMBER OF SHOTS = 5000 



* * 

* QCDINS version 2.0 * 

* * 



NEW 


MAXIMUM 


WEIGHT = 


2 


.4520378923431534E- 


-05 


NEW 


MAXIMUM 


WEIGHT = 


6 


.0858601240747906E- 


-02 


NEW 


MAXIMUM 


WEIGHT = 





.1289870599803515 




NEW 


MAXIMUM 


WEIGHT = 





.3662942637479250 




NEW 


MAXIMUM 


WEIGHT = 





.5705336816482935 




NEW 


MAXIMUM 


WEIGHT = 





.8061467165013312 





INITIAL SEARCH FINISHED 



OUTPUT ON ELEMENTARY PROCESS 



NUMBER OF EVENTS 
NUMBER OF WEIGHTS 
MEAN VALUE OF WGT 
RMS SPREAD IN WGT 
ACTUAL MAX WEIGHT 
ASSUMED MAX WEIGHT 





5000 
2.8721E-02 
6.4377E-02 
7.3286E-01 
8.0615E-01 



PROCESS CODE IPROC = 
CROSS SECTION (PB) = 
ERROR IN C-S (PB) = 
EFFICIENCY PERCENT = 



17600 
28.72 
0.9104 
3.563 



EVENT 



1: 27.50 GEV/C E+ 



ON 820.00 GEV/C P 



PROCESS: 17600 



SEEDS: 106645412 k 135135176 STATUS: 100 ERROR: WEIGHT: 0.2872E-01 
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—INITIAL STATE— 



IHEP ID 


IDPDG 


1ST 


M01 


MO 2 


DAI 


DA2 


P- 


-X 


P- 


-Y 


P- 


-Z 


ENERGY 


MASS 


1 E+ 


-11 


101 


o 





4 








00 


Q 


00 


27 


50 


27 


.50 


00 


2 P 


2212 


102 

















.00 





.00 


-820. 


00 


820 


.00 


0.94 


3 CMF 





103 


1 


2 











.00 





.00 


-792. 


50 


847 


.50 


300.33 








— 


—INITIAL STATE 




















IHEP ID 


IDPDG 


1ST 


M01 


M02 


DAI 


DA2 


P- 


-X 


P- 


-Y 


P- 


-Z 


ENERGY 


MASS 


4 E+ 


-11 


101 


1 














.00 





.00 


27. 


.50 


27 


.50 


0.00 











-HARD SUBPROCESS- 


— 
















IHEP ID 


IDPDG 


1ST 


M01 


M02 


DAI 


DA2 


P- 


-X 


P- 


-Y 


P- 


-Z 


ENERGY 


MASS 


5 E+ 


-11 


121 


7 


9 


21 


9 





.00 





.00 


27. 


.50 


27 


.50 


0.00 


6 GLUON 


21 


122 


7 


20 


22 


19 





.00 





.00 


-50. 


37 


50 


.38 


0.75 


7 HARD 





120 


5 


6 


9 


20 





.53 





.99 


-22. 


98 


77 


.98 


74.51 


8 INST 





3 


7 














.23 


9 


.33 


-43. 


54 


47 


.01 


15.07 


9 E+ 


-11 


123 


7 


5 


26 


5 


8 


.20 


-7 


.47 


3. 


.48 


11 


.63 


0.00 


10 UBAR 


-2 


124 


7 


14 


27 


17 


-8 


.43 


-1. 


,86 


17. 


.19 


19 


.24 


0.31 


11 DQRK 


1 


124 


7 


12 


31 


13 


-0 


.07 





.52 


-7. 


.65 


7 


.67 


0.32 


12 GLUON 


21 


124 


7 


13 


33 


11 





.16 


-0 


.28 


-0, 


.64 


1 


.04 


0.75 


13 DBAR 


-1 


124 


7 


11 


35 


12 


-0 


.67 


-0 


,12 


-2. 


65 


2 


.75 


0.32 


14 UQRK 


2 


124 


7 


15 


37 


10 





.01 


4 


.97 


-11. 


73 


12 


.75 


0.32 


15 GLUON 


21 


124 


7 


16 


39 


14 





.53 





,62 


-2. 


.17 


2 


.44 


0.75 


16 GLUON 


21 


124 


7 


17 


41 


15 





.43 


-0 


.43 


-8, 


.44 


8 


.49 


0.75 


17 GLUON 


21 


124 


7 


10 


43 


16 


-0 


.72 





.09 


-2, 


.40 


2 


.62 


0.75 


18 SQRK 


3 


124 


7 


19 


45 


20 


-0 


.11 


2 


.50 


-4. 


.93 


5 


.55 


0.50 


19 GLUON 


21 


124 


7 


6 


47 


18 





.81 





.93 


-1. 


.46 


2 


.06 


0.75 


20 SBAR 


-3 


124 


7 


18 


51 


6 


-0 


.14 





.52 


-1. 


47 


1 


.64 


0.50 


21 ZO/GAMA* 


23 


3 


5 


7 








-8 


.20 


7 


.47 


24. 


02 


15 


.87 


-21.17 



CHECK OF ENERGY-MOMENTUM CONSERVATION IN THE EVENT: 

Sum_i(P_i) [GeV] = ( 4.87E-13 -5.49E-13 1.59E-12 -4.20E-12 ) 
with 60 stable particles in final state contributing. 
OUTPUT ON ELEMENTARY PROCESS 



NUMBER OF EVENTS 
NUMBER OF WEIGHTS 
MEAN VALUE OF WGT 
RMS SPREAD IN WGT 
ACTUAL MAX WEIGHT 
ASSUMED MAX WEIGHT 



1000 
27258 
2.8837E-02 
6.3371E-02 
7.0199E-01 
8.0615E-01 



PROCESS CODE IPROC 
CROSS SECTION (PB) 
ERROR IN C-S (PB) 
EFFICIENCY PERCENT 



17600 
28.84 
. 3838 
3.577 



49 



References 



[1] A. Belavin, A. Polyakov, A. Schwarz and Yu. Tyupkin, Phys. Lett. B 59 (1975) 
85. 

[2] G. 't Hooft, Phys. Rev. Lett. 37 (1976) 8; Phys. Rev. D 14 (1976) 3432; Phys. 
Rev. D 18 (1978) 2199 (Erratum); Phys. Rep. 142 (1986) 357. 

[3] S. Moch, A. Ringwald and F. Schrempp, Nucl. Phys. B 507 (1997) 134. 

[4] A. Ringwald and F. Schrempp, Phys. Lett. B 438 (1998) 217. 



[5] A. Ringwald and F. Schrempp, hep-ph/9411217 , in: Quarks '94, eds. D. 
Grigoriev, V. Matveev, V. Rubakov, D. Son and A. Tavkhelidze (World Scientific, 
Singapore, 1995) p. 170. 



[6] M. Gibbs, A. Ringwald and F. Schrempp, hep-ph /9506392 , in: Workshop on 
Deep Inelastic Scattering and QCD (DIS 95), eds. J.-F. Laporte and Y. Sirois 
(Ecole Polytechnique, Paris, 1995) p. 341. 

[7] T. Carli, J. Gerigk, A. Ringwald and F. Schrempp, |hep-ph/9906441| , in: Monte 
Carlo Generators for HERA Physics, eds. A.T. Doyle, G. Grindhammer, G. 
Ingelman and H. Jung, DESY-PROC-1999-02 (DESY, Hamburg, 1999) p. 329. 

[8] G. Marchesini et al, Comput. Phys. Commun. 67 (1992) 465. 

[9] S. Moch, A. Ringwald and F. Schrempp, |hep-ph/970640C , in: Deep Inelastic 
Scattering and QCD (DIS 97), eds. J. Repond and D. Krakauer (American 
Institute of Physics, Woodbury, New York, 1997) p. 1007. 

[10] A. Ringwald and F. Schrempp, |hep-ph / 9706399| , in: Deep Inelastic Scattering 
and QCD (DIS 97), eds. J. Repond and D. Krakauer (American Institute of 
Physics, Woodbury, New York, 1997) p. 781. 

[11] S. Aid et al., HI Collaboration, Nucl. Phys. B 480 (1996) 3. 

[12] S. Aid et al., HI Collaboration, Z. Phys. C 72 (1996) 573. 

[13] T. Carli and M. Kuhlen, Nucl. Phys. B 511 (1998) 85. 

[14] S. Moch, A. Ringwald and F. Schrempp: in preparation 

[15] A. Ringwald and F. Schrempp, Phys. Lett. B 459 (1999) 249. 

[16] C. Caso et al., Particle Data Group, Eur. Phys. J. C 3 (1998) 1. 

[17] H. Jung, private communication. 

[18] T. Sjostrand, Comput. Phys. Commun. 82 (1994) 74. 

[19] V.V. Khoze and A. Ringwald, Phys. Lett. B 259 (1991) 106. 

[20] J. Verbaarschot, Nucl. Phys. B 362 (1991) 33. 



50 



[21] G. Grindhammer, private communication. 
[22] A. Mueller, Phys. Rev. D 2 (1970) 2963. 

[23] M. Gibbs and B. Webber, Comput. Phys. Commun. 90 (1995) 369. 

[24] T. Carli, private communication. 

[25] J. Owens, Phys. Lett. B 266 (1991) 126. 

[26] R. Kleiss and W. Stirling, Nucl. Phys. B 385 (1992) 413. 

[27] M. Seymour, private communication. 

[28] J. Bromley et al., in: Future physics at HERA, eds. G. Ingelman, A. De Roeck 
and R. Klanner (DESY, Hamburg, 1996) Vol. 1, p. 611. 



51 



